Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.
N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Ce rapport devra être mis à disposition et partagé avec nous sous sa forme compilée (html sous forme de Github Pages ou à défautPDF) dans votre dépôt github.
Les parties 1 à 4 seront notées pour l’évaluation du module 4, les parties 5 et 6 pour le module 5.
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
mkdir ~/EVALUATION
cd ~/EVALUATION
SRR10390685 grâce à l’outil sra-tools [1]module load sra-tools
fasterq-dump --version
# "fasterq-dump" version 2.10.3
srun --cpus-per-task 8 fasterq-dump -S -p SRR10390685 --outdir . --threads 8
gzip SRR10390685_1.fastq
gzip SRR10390685_2.fastq
zcat SRR10390685_1.fastq.gz | echo $((`wc -l`/4))
zcat SRR10390685_2.fastq.gz | echo $((`wc -l`/4))
ou
module load seqkit
seqkit stat SRR10390685_*.fastq.gz
# file format type num_seqs sum_len min_len avg_len max_len
# SRR10390685_1.fastq.gz FASTQ DNA 7,066,055 1,056,334,498 35 149.5 151
# SRR10390685_2.fastq.gz FASTQ DNA 7,066,055 1,062,807,718 130 150.4 151
Les fichiers FASTQ contiennent 7 066 055 reads.
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
gzip -d GCF_000009045.1_ASM904v1_genomic.fna.gz
module load seqkit
seqkit version
# seqkit v0.14.0
zcat GCF_000009045.1_ASM904v1_genomic.fna |seqkit stat
# file format type num_seqs sum_len min_len avg_len max_len
# GCF_000009045.1_ASM904v1_genomic.fna FASTA DNA 1 4,215,606 4,215,606 4,215,606 4,215,606
La taille de ce génome est de 4 215 606 paires de bases.
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
zgrep -v "^#" GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '($3 == "gene")' |wc -l
# 4448
4 448 gènes sont recensés dans le fichier d’annotation.
module load fastqc
fastqc --version
# FastQC v0.11.9
mkdir QC
srun --cpus-per-task 8 fastqc SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc SRR10390685_2.fastq.gz -o QC/ -t 8
module load multiqc
multiqc --version
# multiqc, version 1.9
srun multiqc -d QC -o QC
La qualité des bases me paraît … car … comme le montre …
Lien vers le rapport MulitQC
Oui/Non car …
(R1+R2) / Taille du génome
La profondeur de séquençage est de : 502.6898187 X.
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp qui vous semblent adéquats et justifiez-les.
module load fastp
fastp --version
# fastp 0.20.0
mkdir FASTP
srun --cpus-per-task 8 fastp --in1 SRR10390685_1.fastq.gz --in2 SRR10390685_2.fastq.gz --out1 FASTP/SRR10390685_1.fastq.gz --out2 FASTP/SRR10390685_2.fastq.gz --html FASTP/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json FASTP/fastp.json
seqkit stat FASTP/SRR10390685_[12].fastq.gz.
Les paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| –cut_mean_quality | 30 | pour un score moyen dans la fenêtre glissante > 30 |
| –cut_window_size | 8 | pour une taille de fenêtre glissante de 8 |
| –length_required | 100 | pour ne garder que les reads de taille > 100 |
| –cut_tail | pour faire partir la fenêtre de l’extrémité 3’ du read |
Ces paramètres ont permis de conserver 6 777 048 reads pairés, soit une perte de 4.0900757 % des reads bruts.
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [3] et samtools [4].
module load samtools
samtools --version
# samtools 1.10
# Using htslib 1.10.2
module load bwa
bwa
# Version: 0.7.17-r1188
srun bwa index GCF_000009045.1_ASM904v1_genomic.fna.gz
mkdir MAPPING
srun --cpus-per-task=4 bwa mem GCF_000009045.1_ASM904v1_genomic.fna FASTP/SRR10390685_1.fastq.gz FASTP/SRR10390685_2.fastq.gz -t 3 | samtools view -hbS - > MAPPING/SRR10390685.bam
srun samtools flagstat MAPPING/SRR10390685.bam > MAPPING/SRR10390685.bam.flagstat
srun samtools sort MAPPING/SRR10390685.bam -o MAPPING/SRR10390685_sorted.bam
srun samtools index MAPPING/SRR10390685_sorted.bam
samtools view -f 4 -c MAPPING/SRR10390685.bam
# 744540
744 540 reads ne sont pas mappés.
Le gène trmNF est indispensable pour la suite des analyses. Pour vérifier s’il a bien été séquencé, vérifiez s’il est couvert et à quelle profondeur grâce à l’outil bedtools [5]:
module load bedtools
bedtools --version
# bedtools v2.29.2
grep trmNF GCF_000009045.1_ASM904v1_genomic.gff | awk '$3=="gene"' > trmNF.gff3
bedtools intersect -a MAPPING/SRR10390685_sorted.bam -b trmNF.gff3.gz -f 0.5 > SRR10390685_on_trmNF.bam
samtools view -c SRR10390685_on_trmNF.bam
2 801 reads chevauchent le gène d’intérêt.
1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.
4. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.
5. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France